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Introduction 

The evolution of high-speed computers and so- 
phisticated display devices has encouraged the devel- 
opment of advanced algorithms for manipulating and 
displaying multidimensional data. In particular, the 
area of computer-aided geometric modeling (CAGM) 
has advanced significantly in recent years. In CAGM, 
computational geometry and computer graphics are 
combined to give mathematical and graphical repre- 
sentations of curves, surfaces, and volumes. A wide 
variety of applications may be found in the mathe- 
matical representation of physical phenomena, such 
as meteorological data, and in the design of aircraft. 

Often it is desirable to smooth three-dimensional 
surface data consisting of two independent variables 
( x and y) and one dependent variable which contains 
random noise because of errors in measurement or 
calibration. The purpose of smoothing the surface is 
to obtain statistically representative values of the de- 
pendent variable. One typical approach to smoothing 
surfaces with splines is to first smooth along rows of 
data and then smooth along columns of data (ref. 1). 
If, however, the data are not aligned in rows and 
columns, then additional procedures must be applied, 
such as the use of a triangularization method to in- 
terpolate to a rectangular grid before smoothing. A 
difficulty that can arise with either method is that 
changes in trends in the data can induce undesirable 
oscillations in the smoothed spline surface. 

One way to reduce or eliminate the oscillations in 
the spline surface is to apply tension to the surface. 
Applying mathematical tension to a spline surface 
is analogous to grasping the opposite edges of a 
membrane and stretching the membrane to remove 
wrinkles. In reference 2, Spath developed the rational 
spline for both curve and surface interpolation. The 
rational spline is a cubic function which has tension 
parameters in the denominators of the cubic terms. 
More recently, Frost and Kinzel (ref. 3) developed 
an algorithm for automatically adjusting the tension 
parameters in curve-interpolating rational splines. 
The Frost and Kinzel method was combined with 
constrained least squares by Schiess and Kerr (ref. 4) 
to give an algorithm for rational-spline smoothing 
of curves. Finally, Schiess (ref. 5) developed two 
algorithms for rational-spline interpolation of surface 
data given on a rectangular grid. 

The present paper presents an algorithm for 
smoothing surface data with bivariate rational 
splines having multiple tension parameters. The 
multiple tension parameters allow for local con- 
trol of tension on the smoothing surface. Equa- 
tions are derived to ensure continuity of the deriva- 
tives at the knots. Smoothing is accomplished by 


finding the weighted least-squares estimates of the 
bivariate rational-spline coefficients. The capabili- 
ties of the rational-spline smoothing algorithm are 
demonstrated on terrain elevation data. 
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4 by 4 matrix of coefficients 

coefficients of multivariate rational 
spline in the subregion Rij (fc, l — 

1,2, 3, 4) 

coefficients in equation for FX^j 

coefficients in equation for FY % j 

coefficients in univariate rational 
spline on interval i ( k = 1, 2, 3, 4) 

2 MN by 4 MN matrix of constraint 
cofactors 

difference between consecutive values 
of ^-coordinate of knot, \ — X{ 

difference between consecutive values 
of ^coordinate of knot, Vj+i — Vj 

m by 4MN matrix of rational-spline 
cofactors 

function value at the data point 
(x r , 2/r ) 

function value at (x 2 -, yj) 

partial derivative of /^(x, y) with 
respect to x evaluated at (x z -, yj) 

partial derivative of fij{x,y) with 
respect to y evaluated at (x^,yy) 

partial derivative of /^-(x,y) with 
respect to x and y evaluated at (x t *,yy) 

rational spline on the subregion R^j 

4 by 4 matrix of functions gik{x) 

functions of x used in rational-spline 
representation on interval i ( k — 

1 , 2 , 3 , 4 ) 

4 by 4 matrix of functions hji(y) 

functions of y used in rational-spline 
representation on interval j (l = 

1,2, 3, 4) 

2MiV-element column vector of 
Lagrange multipliers 

number of knots along x-axis 


m number of data points 

N number of knots along y - axis 

Pi tension factor for interval i 

q, j tension factor for interval j 

R tJ rectangular subregion in xz/-plane 

defined by X{ < x < x t -+i and 
Vj <V< Vj+i 

S ^ 4 by 4 matrix of function and deriva- 

tive values 

s ij 16-element column vector of unknowns 

r, s, t , u variables used to define rational spline 

Vj(x) 4-element column vector of functions 

of x 

vw ij 16-element vector of functions of x and 

y 

V estimated variance of measurement 

error 

W m by m diagonal weighting matrix 

w j(y) 4-element column vector of functions 

of y 

x, y independent variables 

Z ra-element column vector of function 

values F r 

Z 4AL/V-element column vector of un- 

known parameters 

z dependent variable in univariate 

rational spline 

A x , Ay differences with respect to x and y 
variables 

Subscripts: 

f, j quantity at the knot (x^,?/y) 

k , l general indices 

r data point index 

Superscripts: 

T matrix transpose 

— 1 matrix inverse 

A prime indicates first derivative with respect to 
the independent variable. A double prime indicates 
second derivative with respect to the independent 
variable. A bar over a quantity denotes that quantity 
at a knot. 


Problem Statement 

Let the three-dimensional data (x r , 2/r, F r ) be 
given, where r = 1,2, The variables x and 

y are the independent variables. The paired val- 
ues (x r , y r ) lie in a finite, bounded region of the 
xy- plane but do not need to be uniformly scattered 
in the region or lying on grid lines. The values F r of 
the dependent variable represent values of a function 
measured at the points (x r , y r )\ the measured val- 
ues are assumed to be corrupted by random error of 
unknown statistical characteristics. 

The objective is to fit, in the weighted least- 
squares sense, a bivariate rational spline to the data. 
This requires the selection of M knot locations along 
the x-axis and N knot locations along the y ~ axis so 
that xi < x 2 < ... < xm and y\ < y 2 < ••• < Vn- 
The knot locations need not be equally spaced. All 
the data must lie in the region defined by {xi < x < 
XAf , y i < y < y^}, and several data points should 
lie in each subregion R^j = {x\ < x < x^ +1 , yj < 

y < yj+ 1 ). 

Rational Splines 

In this section both the univariate and bivariate 
rational splines are described. Because the bivariate 
rational spline is the tensor product of two univariate 
rational splines, the characteristics of the univariate 
rational spline are discussed first. 

Univariate Rational Spline 

Let xi be the abscissas of the knots of the uni- 
variate rational spline, where i = 1,2, ...,M and 
x\ < X 2 < ... < xm, and let z be the value of 
the spline at x. The rational spline on interval 
i {i = 1, 2, ..., M — 1) is defined in references 2 and 3 
to be 

4 

2 = 52 9ik{ x ) i x i < x < x i+ 1) (1) 

k= 1 

where are unknown coefficients, 

u 3 

9i lOO = « 9ii{x) = — — 

Pit + 1 

t 3 

9i2( x ) = t gi 4 {x) = — 

PiU + 1 

where 

X^_j_ i x 

* =j v 


2 



dxi 
dxi = 

and Pi is the tension parameter for interval i 

Equation (1) is defined for all values of the inde- 
pendent variable x in the data range if the tension 
parameter pi is restricted to Pi > —1. If Pi is set to 
zero, equation (1) reduces to a cubic-spline function. 
As pi increases from zero, the cubic terms decrease 
in magnitude and the function tends to the equation 
of the line joining the knots at x z * and x^i. Because 
a distinct, independent tension parameter is associ- 
ated with each interval, the behavior of the function 
in each interval may be locally controlled. 

Evaluation of equation (1) for each subinterval re- 
quires knowledge of the four coefficients c^, 
and c^. Thus, for M data points (equivalently, for 
M — 1 subintervals), 4 M — 4 coefficients must be de- 
termined. Spath (ref. 2) reduces the magnitude of 
this problem by writing the coefficients in terms of 
the values of the function and its first derivative at 
the knots. End conditions are applied to the first 
derivative, and equations ensuring the continuity of 
the second derivative at the interior knots are de- 
rived. For the interpolation problem, this derivation 
yields a system of M — 2 equations for the M — 2 
unknown interior first derivatives. Frost and Kinzel 
(ref. 3) extend Spath’s approach by allowing for three 
different end conditions and by developing an itera- 
tive method for determining the tension parameters. 
Tension parameters are found so that the interpo- 
lating rational spline deviates from the line joining 
knots by a prescribed value. 

Another approach to determining a smooth fit 
to the data has been presented by Schiess and Kerr 
(ref. 4) in deriving a least-squares univariate rational- 
spline approximation. The rational spline is reformu- 
lated in terms of the unknown spline function and 
its second derivative at the knots. Smoothness is 
ensured by imposing the constraints that the first 
derivatives are continuous at the interior knots. This 
approach leads to a constrained least-squares prob- 
lem in the 2 M values of the unknown function and 
its second derivative at the knots. 

Bivariate Rational Spline 

Let a given set of M by N knot points in three 
dimensions be represented by (xj, t/y, F t y), where 
i = 1,2,...,M and j = 1,2,..., iV. The independent 
variables are assumed to be ordered (x\ < X 2 < • •• < 
xm and yi < 2/2 < ••• < Vn) and form a rectangular 
grid, but are not necessarily equally spaced. 


The multiple-tension-parameter bivariate rational 
spline on the subregion R % j defined by Xi < x < 
*i + 1 (i = 1 , 2 ,..., M - 1) and yj < y < j/y +1 (j = 
1,2, ..., N — 1) is defined in reference 2 by 

4 4 

fij(x>V) = E H a ijkl 9ik( x ) h jl(v) (2) 

k=l 1=1 

where gik(x) and Pi are the same as for the univariate 
spline, aij^i are unknown coefficients, 

s 3 

hji(y) = s h jz {y) = —-^ 

r 3 

M») = r 

where 

. yy+i - y 
dyj 

V-Vj . 

r = — - — - = 1 — s 

d yj 

dyj = y j+1 - yj 

and qj is the tension parameter for interval j. 

As defined by equation (2), the bivariate rational 
spline on each subregion Rij is a function of 2 tension 
parameters ( Pi and qj) and 16 coefficients (a^/). 
The coefficients are to be determined so that the 
rational spline and its first and second derivatives are 
continuous over the entire region. The M + iV-2 
tension parameters may be adjusted individually; 
each parameter affects the behavior of the rational 
spline in a strip parallel to either the x-axis (for qj) 
or the y - axis (for pi). 

Since 16 coefficients are needed on each subregion, 
a total of 16 ( M - l)(Af — 1) coefficients must be 
determined to define the entire rational spline. For 
example, for a 30 by 30 grid (M = N = 30), a total of 
13 456 coefficients are needed. Spath (ref. 2) reduces 
the actual number of unknown quantities by writing 
the coefficients as linear combinations of the values 
of the function and its derivatives at the grid points. 
A similar approach is taken in this paper. 

Let FXij and FY ij be the first derivatives of 
fij{x,y) with respect to x and respectively, and 
FXYij be the cross derivatives, all evaluated at the 
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point (xj,j/y). For the subregion Rjj , define the 4 by 
4 matrices 


' Fii 

FYii 

F 

FY iu +» ■ 

FXij 

TxYij 

F* <(i+D 

FXY iU+i) 

F(i+l)j 

FYi i+ib 


FY (i+i)u+i) 

-FX {i + 1) 3 

FXY {iJri)j 


FXY ( i + 1 )( J + 1 ) . 


A ti 


“ijli 
a ij 21 
a ij31 


a tjl2 

a ij22 
a ij 32 


L Q, 2 j41 a ij42 


a ij 13 
a ij23 
a ij33 
a ij 43 


a ij 14 
a ij 24 
a ij 34 
a ij44 - 


G, = 


9.3^) 9i4(*i) 


^3(^1) ^4(^1) 


9.1 (*i) 9.2 (*i) 

9il( 5 «) 9' 2 ( s i) »»4' 

9.1 (^.+ 1 ) 9i2( J i + l) 9i3(2.+ l) 9,4( i i+l) 
L9- 1 (*«+l) 9' 2 (*,+i) 9' 3 (ii+ 1 ) 9- 4 (i.+ l) 


H i = 


- 1 0 1 

1 1 3 + p t 

dx± dx, dxi 

0 1 0 

\ l n 

dx± dx± 

0 ■ 
0 
1 

3 + Pl 

dx± J 


h jl{yj) 

fcj2 iVj) 

3^ 

CO 

1 h ]4 (yj) 

h 'jM) 

h 'j 2^j) 

h j3^ y J J 

1 h j 4 ( %) 

h j2(yj+ 1) 

h j3(Vj+l) h j4(Vj+l) 

+ h 'j2^3+ 1) 

h j3(yj+i) ^ 4 (yj+i)- 


" 1 

0 

1 

0 ■ 

_ 1 

i 

3 + Qj 

“3 7T 

0 

0 

1 

0 

1 

1 

L 

1 

3y“ 

0 

3 + qj 
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Thus, in matrix notation, 

S tj = GjAjj-HJ (3) 

Equation (3) can be verified by differentiating equa- 
tion (2), as appropriate, and evaluating the results at 
the corners of the subregion. Note that Gj and H j 
depend on the grid spacing and tension parameters 
but not on the function values. 


Since the matrices Gj and Hy are nonsingu- 
lar, equation (3) can be solved for the matrix of 
coefficients: 


A ij = G- 1 S ij {Hjr 1 (4) 

Therefore, for any subregion the 16 coefficients can 
be determined from the values of the function, its 
first derivatives with respect to x and t/, and its 
cross derivatives at the four corners of the subregion. 
Therefore, a total of 4 MN function and derivative 
values are needed to calculate the coefficients. 

Rational-Spline Smoothing 

In this section an algorithm for finding the 
surface-smoothing rational spline in terms of the 
function and its derivatives is presented. This algo- 
rithm is for the general case of M— 1 values of the ten- 
sion parameters Pi {i = 1, 2, ..., M — 1) and N —l val- 
ues of the tension parameters qj (j ~ 1, 2, ..., N — 1). 
An algorithm for the single-tension-parameter ra- 
tional spline is not presented because that rational 
spline is a special case of the multiple-parameter 
spline. 

Vector Formulation 

Let the knot locations (ay, yj) and tension param- 
eters Pi and qj (for i = 1,2, ..., M and j — 1, 2, ..., N) 
be given. The knots do not need to be equally spaced, 
but they must form a rectangular grid. Let the data 
(ay-, y r , F r ) for r = l,2,...,m also be given. Least- 
squares estimation of the 4MN unknown function and 
derivative values requires that there be more data 
points than unknowns, or that m > 4 MN. Further, 
there should be several data points in each region 
ityy. __ 

In order to find least-squares estimates of F t j , 
FX{j , FY ij , and FXY it is necessary to reformu- 
late the problem. Using equation (4) in equation (2) 
results in the rational spline being written as 


fij(x,y ) = vf(x) SijWj{y) (5) 


where 


Vi(x) = (G. 1 ) T [&i(x), g i2 (x ), & 3 (x), g l4 {x) } T 

( 6 ) 

and 


w j(y) = (Hy *) T [hji(y), h j2 (y), h j3 {y), h j4 (y)] T 

( 7 ) 
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In terms of individual entries in v^(x), wy(y), and 
S {j, equation (5) can be written as 

4 4 

fij(x,y) = ^2 53 v **( x ) s ijkl Wjdv) 

k= 1/=1 
4 4 

= E£ ▼»*(*) w i/(y) SyW (8) 
fc=n=i 

or 

fij{x,y) = I™,/ Sy 

with the two 16-element column vectors defined by 
vwii = ki(z)wji(y), W f i(*)u>j2(y), Wji(x) 


(y ) » — » v »4 (a?)Wj4 (y)] 3 


and 


s ij — [SijUi Sijl2* 1 3 , Siyi 4 , Sjj* 21 »***> *^zj44 ] 

As presented in equation (8), f l j{x, y ) is a linear func- 
tion of the vector Sjj of unknown function and deriva- 
tive values at the knots and therefore is amenable to 
least-squares estimation. 

One minor modification which is made to simplify 
the use of this formulation is to rearrange the entries 
in s ij (and corresponding entries in vw jj) so that 
all the unknowns corresponding to one knot are in 
consecutive entries. The rearrangemen t is cho sen so 
that the entries are in the order F^-, FX t y, FYjj, 
and FXY ij. 


(with the current notation and some rearrangement). 
Constraining the second derivatives with respect to 
x to be continuous at the interior knots results in the 
(M — 2)iV equations 


CXi_i FX(^i)j 4- [(2 + Pi-\)CXi_i 

+ (2 + Pi )CX l \FX ij + CXiFX^j 
3 + Pi- 1 


dxj^ i 

3 + Pi 


CXi _ i A X F 


dXl cx *• = 0 
(*' = 2, 3, M — 1; j = 1,2,..., AT) 


( 9 ) 


where 


CX = 


pj + 3p z + 3 

[(2 + Pi) 2 - l}dxi 


^x F ij — F^ +1 )y F{j 

Similarly, continuity of the second derivative with 
respect to y at the interior knots yields the additional 
M(N — 2) equations 

CVj-i FY i{j _ D + [ (2 + qj-JCYj-i 
+ (2 + q^CYjUFYij- + CYjFY W) 

3 + Qj-i 


dyj - 1 


CYj— l 


Constraint Equations 

In the formulation used herein, the continuity of 
the first derivatives and cross derivatives at the in- 
terior knots is ensured because these derivatives are 
parameters to be estimated in the rational-spline for- 
mulation. It is also desirable that the second deriva- 
tives be continuous at the interior knots. Further- 
more, for the knots on the boundaries, natural spline 
conditions can be applied (i.e., the second derivatives 
are zero on the boundary). Both of these conditions 
can be added to the least-squares formulation as con- 
straint equations which are linear in the unknown 
quantities. 

For the continuity of the second derivatives at the 
interior knots, the equations derived in reference 5 
are used. Equations (9) and (10) of this report are 
the same as equations (5) and (6) in reference 5 


3 -f q, — 

- — ^ CYj A y Fjj = 0 
dyj 3 y 13 

(* = 1,2,...,M; j = 2, 3, ..., N — l) (10) 

where 

q) + Mj + 3 

J [(2 + q j ) 2 -l]dy j 

^y F ij = F i(j+ 1 ) — F ij 

The natural spline boundary conditions require 
that the second derivatives with respect to x be zero 
at the boundaries, conditions which are equivalent to 
requiring the rational spline be linear on and outside 
the boundaries. These conditions result in two sets 
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of constraint equations. Along the boundary x = xi, 
the equations are 

3 -zr r * - + < 2 + 

+ FX 2 j = ^ U= 1 , 2 , ( 11 ) 


Along the boundary x = the equations 


3 + PM-ln 


3 + Pm-It; 

— j t Mj 

ax M _i 


+ + (2 + pm-i) FXjtfj — 0 

ti = 1,2,..., AT) 


(12) 


On the boundaries y = y\ and y — yj y, the second 
derivative with respect to y is zero. These conditions 
lead to two additional sets of constraint equations. 
For the boundary y = i/i, equation (13) applies: 


^-^•1 - ^-^2 + (2 + < 71)^1 
dm dy i 

+ FF t2 =0 (* = 1, 2, Af) (13) 

Along the boundary y = y N , the constraint equations 
are 


that Z = (Fi, F 2 , Frn)_ T . Defi ne the m by 
4MN matrix E of cofactors F ZJ , FX^-, FY ij , and 
FXYij (* = 1,2, ,..,M and j = 1,2, ...,7V) as given 
in equation (8). F rom the defini tion of Z, the quadru- 
plets (F t *y, FXij, FY ij , FXy^) are ordered so that 
the subscript z changes faster than the subscript j. 
Equivalently, the values in Z are ordered so that all 
the values in region Rn are first, those in region F 21 
are second, and so forth; the ordering within regions 
is unimportant. With this ordering, E has the gen- 
eral block structure 



E = 


L 




where only the entries in blocks are nonzeroes. Each 
block is eight columns wide. With these defini- 
tions, equation (8) evaluated at (x r , y r , F r ) (r — 
1, 2, ..., m) can be written 


3 + 97V— l 
dVN - 1 


F i(jV-l) 


3 


+ QN-l 

dyN - 1 


FiN 


+ Fy^jv-i) 


+ l2 + q N -i) FY iN = 0 (* — 1,2, ..., M) (14) 


The (M — 2)7V + M(N — 2) constraints at the 
interior knots and 2 M + 2N boundary constraints 
yield a total of 2 MN constraint equations. The 
derivation of the constraint equations resulting from 
the boundary conditions is given in the appendix. 
In the next section these equations are incorporated 
into the least-squares fit of the rational spline to the 
data. 


Constrained Weighted Least Squares 

The rational-spline surface approximation is 
found through solution of a constrained weighted 
least-squares problem. The solution to this prob- 
lem is most easily seen with the following matrices 
defined. Let Z be the 4M/V-element column vector 
of unknowns such that Z = (Fn, FXq, FY n, 
FXYn, * 21 , FX 2h ...,FXY MN ) T . Let Z be 
the m by 1 vector of known function values F r 
at the data points (x r , y r ) (r = 1,2, ...,m) such 


Z = EZ (15) 

Let D be the 2MNby 4M/V matrix of cofactors of 
the constraint equations (eqs. (9) to (14)). Thus, the 
constraint equations can be written in matrix form 
as 

DZ = 0 (16) 

Each row of D has only four (corresponding to 
eqs. (11) to (14)) or six (eqs. (9) and (10)) nonzero 
entries. 

Finally, let W be an m by m diagonal weight- 
ing matrix containing the weights on the diagonal. 
Then, the constrained weighted least-squares prob- 
lem requires minimization of (Z - EZ) T W(Z — EZ) 
with respect to Z such that DZ = 0. With L defined 
as a 2AFV-element column vector of Lagrange multi- 
pliers, the constrained problem can be rewritten as 
the unconstrained problem, that js, minimization of 
(Z - EZ) r W(Z - EZ) + 2L r DZ with respect to Z 
and L. The solutions of this minimization problem 
(ref. 6) are 

Z = (E r WE) _1 E r WZ - (E r WE)“ 1 D t L (17) 
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and 

L = [D(E t WE)' ! D t ]- 1 D(E t WE)' 1 e t wz 

(18) 

The solutions in equations (17) and (18) exist and 
are unique if the matrix inverses in the equations 
exist. Furthermore, a sufficient condition such that 
the solutions in these two equations yield_a unique 
global minimum of (Z — EZ)-^W(Z — EZ) subject 
to the constraint DZ — 0 is Z T E T WEZ > 0 for all 
Z > 0 with DZ = 0 (ref. 6). 

For the user interested in an overall measure of the 
error of the fit to the data, the estimated variance 
of the measurement error is easily calculated. The 
estimated variance V is defined by (ref. 6) 

V = m-iMN (Z-EZ) T W(Z-EZ) (19) 

Selection of Tension Parameters 

Although an automated procedure for adjusting 
tension parameters is available for univariate rational 
splines (refs. 3 and 4), no such procedure has been 
devised for bivariate rational splines. Instead, it is 
recommended that the user visually examine three- 
dimensional or contour plots of the data and the 
rational-spline smoothed surface and then use engi- 
neering judgment to select tension-parameter values. 
This approach is outlined here. 

First, the user should obtain a plot of the orig- 
inal data. This plot will show both general trends 
and local anomalies of the data. Second, the user 
should compute the smoothing bicubic-spline surface 
by calculating the smoothing rational spline with all 
the tension parameters set to zero. Comparison of 
the data and the bicubic-spline plots will indicate re- 
gions where the cubic spline exhibits undesirable or 
exaggerated hills and valleys. Using small to mod- 
erate tension-parameter values (say, from 1 to 10) 
for those regions, the user can recalculate a smooth 
rational spline. In this way, after a few trials with 
different tension values, a rational-spline surface can 
be obtained which the user considers representative 
of the data. 

Implementation Considerations 

In implementing the rational-spline surface 
smoother, the number and locations of the knots need 
to be considered. For the given m data points, the 
number of knots chosen (M and N) is restricted by 
the relationship 

AMN < m (20) 


since there are 4 MN unknowns. The knot locations 
need to be selected so that all the data lie within 
the region {xi < x < x m , y\ < y < y M } or 
on the boundaries. A few data points must lie in 
each subregion Since there are ( M — 1)(JV — 1) 
subregions, dividing each side of equation (20) by 
( M — 1 )(N — 1) indicates that there must be more 
than 4 MN/(M — 1) (TV — 1) data points either in each 
subregion or on the boundaries of the subregion. 
Note that points on a boundary between two subre- 
gions count towards the total number of data points 
for each of the two subregions. In the terrain eleva- 
tion example given later, m — 121 and M — N — 5. 
Thus, equation (20) becomes 100 <121 and there 
must be more than (4) (5) (5)/ (4) (4) = 6.25 data 
points per subregion; that is, each subregion must 
contain at least 7 data points. As constructed, there 
are 9 data points within or on the boundaries of each 
subregion. 

The two matrices to be inverted in the solution 
of equations (17) and (18) can be large — E T WE is 
AMNby 4MAand D(E T WE)” 1 D T is 2MNby 2 MN. 
For example, for M — N = 5 these matrices are 100 
by 100 and 50 by 50, respectively. Calculation of 
such large inverses with a general matrix inversion 
method, such as that of Gauss- Jordan, may lead 
to loss of precision. However, both matrices are 
symmetric and positive definite. Therefore it is 
highly recommended that a Cholesky method (ref. 7) 
be used since it takes advantage of the characteristics 
of the matrices. 

Example 

The data chosen to illustrate the rational-spline 
surface smoother consist of terrain elevation mea- 
surements on a square area 1000 ft on each side. 
The measurements were taken at 100-ft intervals in 
both the x- and the ^direction. This yields a total 
of m = 121 measurements. Elevations are measured 
to an accuracy of 0.1 ft. 

Figure 1 shows a surface plot of the terrain 
elevation data. The elevations range from 12.0 ft 
(at x — 300 ft, y = 300 ft) to 31.0 ft (at x = 700 ft, 
y — 500 ft). As shown in figure 1, the original data 
are relatively smooth. For the purposes of this ex- 
ample, the original elevation measurements are per- 
turbed by addition of normally distributed error hav- 
ing a mean of zero and a standard deviation of 1 ft. A 
surface plot of the resulting measurements is shown 
in figure 2. It is to these measurements that the sur- 
face smoother is applied. 

Knots are chosen to be located at the points 
0 ft, 250 ft, 500 ft, 750 ft, and 1000 ft along each 
axis; this choice gives a 5 by 5 (N = M = 5) grid of 
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knot locations. As a result, 4 MN = 100 unknown 
function and derivative values must be estimated. 
In the case of zero-tension values (pi = 0 for i = 
1,2, 3, 4 and qj = 0 for j = 1,2, 3, 4), the bivariate 
rational-spline smoother reduces to a bicubic-spline 
smoother. Figure 3 shows a surface plot of the 
bicubic-spline smoothed values of the measurements 
including error. The values plotted in this figure are 
generated through evaluation of the bicubic spline 
at 50-ft intervals along each axis in order to show 
the behavior of the surface at intermediate points. 
Although the bicubic-spline surface is smooth, it 
contains oscillations which are not present in the 
original data. It is these artificial oscillations which 
can be removed or reduced by application of tension. 

Figure 4 shows the rational-spline surface when all 
8 tension parameters have a value of 10. Comparison 
with figure 3 shows that the nonzero tension gives 
the surface a less “rounded” appearance because the 
magnitude of the cubic terms is smaller when the 
tension has a value of 10 than when it has a value of 
0. This tendency for the rational spline to be more 
linear for a tension value of 10 is clearly shown along 
the y = 0 edge for large x. This edge in figure 4 
also has shallower undulations than the same edge in 
figure 3. This reduction in undulations is the purpose 
for applying tension. 

As an example of the application of tension to 
subregions of a surface, only the tension parameters 
P 4 and are set to a value of 10 and all the others are 
set to 0. Figure 5 shows the surface resulting from 
these tension parameters. Comparison of figure 5 
with figures 3 and 4 shows that the surface for the two 
regions (750 ft < x < 1000 ft, 0 ft < y < 1000 ft) and 
(0 ft < x < 1000 ft, 750 ft < y < 1000 ft) in figure 5 
closely resembles the surface for the same regions in 
figure 4. In contrast the remaining surface of figure 5 
more closely resembles the corresponding surface in 
figure 3. Hence, the surface can be smoothed with 


rational splines in which the tension selectively varies 
from region to region. 

Concluding Remarks 

An algorithm for surface smoothing with bivari- 
ate rational-spline functions on a rectangular grid has 
been presented. The rational-spline function com- 
bines the advantages of a cubic function having con- 
tinuous first and second derivatives over the entire 
grid with the advantage of a function having vari- 
able tension. Adjustment of the tension parameters 
in the rational spline allows the user to reduce un- 
wanted oscillations to any desired level across the sur- 
face. The multiple parameters of the rational spline 
provide adjustable parameters for each rectangular 
subregion and thus give the user control over local 
behavior of the surface. 

The terrain elevation example presented illus- 
trates the reduction in undesirable oscillations that 
is possible with a bivariate rational spline which 
smooths the measurements. The example also illus- 
trates that the bivariate rational spline provides the 
capability of controlling local oscillations in the sur- 
face caused by trend variations in the data. 

There is one major disadvantage to using the 
bivariate rational spline. Several computer runs with 
different tension values may be necessary to find 
a smoothing surface satisfactory to the user. It 
is recommended that for the initial run the user 
apply zero tension in order to ascertain regions of 
excessive oscillations. The user can then determine 
appropriate tension values by examining the surface 
plots from a few additional computer runs. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
March 30, 1987 



Appendix 

Derivation of Boundary Constraints 

In this appendix equations (11) to (14), which 
express the constraints along the boundaries, are 
derived. The derivations of all four sets of equations 
are similar. A detailed derivation of equation ( 11 ) is 
given. 

With the definition of the rational spline in equa- 
tion ( 2 ), consider the interval X{ < x < fix 

y, and evaluate fij{x, y) at yj. Thus, hji(yj) = 
hjziVj) = r h j2 (yj) = hjJSj) = 0 , and 


The second derivatives of gik(x ) are 
9<i(*) = s" 2 ( x ) = 0 

„ . , 6«(p t t + l) 2 + 6u 2 pj(p t t + 1) + 2u 3 p? 

I . _ ( Tl = L 


»" 3 (*) = 
= 


dxj(p t t + l ) 3 

6t(Pi u + l) 2 + 6t 2 p^(pjU + 1) + 2t 3 p? 
dx 2 (p^u + l) 3 


(A5) 


Differentiating equation (Al) twice with respect to x, 
applying equations (A5), and evaluating the results 
at Xi and leads to 


fij(x,Vj) = Y2 9ik(x){aijki + a ijki) (Al) 
k= 1 

For simplicity, define bjj^ = + < 2 ^ 3 . Then, 

from equation (Al) and the definitions of ^(x), 


2p 2 + 6 p t + 6 

-epte'W = - J -S 1 


d2 fa 

= 


2p 2 + 6 pj + 6_ 

dx 2 


vij4 


(A 6 ) 


fijiVii Vj ) — 1 T — F z y 

/z'y(^i-hi i Vj ) t>ij2 T ^(i+Dy 


(A2) 


Solve equations (A 2 ) for and 6 ^ 2 , substitute 
into equations (A4), and rearrange the results to 
obtain 


where F{j and are unknown function values 

at the grid points. 

The first derivatives of Qik{x) are 


(2 + Pi)b{j3 b ij4 ~ dx^FXij + F t j F^+^j 

(A7) 


9 f t iW 


1 , ~3u 2 (p^ + 1 ) - u 3 p z 

9% 3 dx^p;* + l ) 2 


J_ / / x 3t 2 (p,M + 1) + t 3 Pi 
d x i dxifau + l ) 2 


(A3) 


Differentiating equation (Al) with respect to x, ap- 
plying equations (A3), and evaluating at X{ and 
yields 


bij3 + (2 + Pi)bij4 — dxiFX^+^j + F {j — F 

(A 8 ) 

Define A x F^j = — Fjj and solve equa- 

tion (A7) for b l j 4 to get 

Wj4 = ~ (2 4- Pi)b{j3 — dx{FX{j + A X F {j (A9) 


Substituting equation (A9) into equation (A 8 ) and 
solving for 6^3 yields 




6 »Jl 

6 »j2 3 + Pil 

d Xi 

dx* dx* 1:7 3 

b ij 1 

b ij2 3 + Pi , 

" *7 + 

dx t + dx, 

FX (*+ib 



(A4) 

where FX z y and are the unknown first 

derivatives with respect to x at interior grid points. 


b ij 3 = 


(3 + Pi)&xF{j (2 -I- p{)dx^FX^j 

(2 + Pi ) 2 — 1 

(A 10 ) 


Substituting equation (A10) into equation (A9) gives 


b ij 4 — 


— (3 + Pi)AiFjj + (2 + Pi)dxjFX( i+l )j + dx^FX^j 
(2 + p t ) 2 - 1 

(All) 
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In order to derive equation (11), evaluate the first 
of equations (A6) at i = 1, substitute for 6^3 from 
equation (A 10), and set the result to zero: 


d 2 fij(xi,yj) _ A 

d* 2 


or 


2p\ 


+ 6pi + 6 

dx ^ 


(3 + p\)A x Fij — dx\FX 2 j - (2 + P\)dx\FX l3 
( 2+ Pl ) 2 - 1 


Substituting for A x Fij and dividing both sides by 
the left factor, which is nonzero, yields 


( 3 + Pi) p 
dx! F 


1 3 ~ 


dxi 


(A12) 


Equation (A12) is identical to equation (11) and 
holds for j = 1, 2, N. 

For the derivation of equation (12), set to zero 
the second derivative with respect to x along the 
boundary x = xj and substitute the second of 
equations (A6) with i = M — 1 and the definition 
of b{j4 in equation (All). 

Equations (13) and (14) are derived analogously 
by differentiation of equation (2) twice with respect 
to y when x = x 1 for fixed i. The result is evaluated 
separately at the boundaries y = y\ and y = y^, 
and appropriate substitutions are made. 
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Figure 2. Terrain elevation data with random error (standard deviation of 1.0 ft) added. 
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Figure 3. Bicubic-spline smoothed terrain elevation surface. 
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